Readout of a qubit coupled to a cavity mode via Jaynes-Cummings interaction

Ondrej Cernotik

23 February 2015

The evolution of a qubit and a thermal cavity coupled via a Jaynes-Cummings interaction is given by the stochastic master equation $$ \mathrm{d}\rho = -i[g(a\sigma_+ e^{i(\omega t+\phi)}+a^\dagger\sigma_- e^{-i(\omega t+\phi)})+\Delta a^\dagger a,\rho]\mathrm{d}t +\kappa(\bar{n}+1)\mathcal{D}[a]\rho\mathrm{d}t+\kappa\bar{n}\mathcal{D}[a^\dagger]\rho\mathrm{d}t +\sqrt{\frac{\kappa}{2\bar{n}+1}}\mathcal{H}[(\bar{n}+1)ae^{-i\delta t}-\bar{n}a^\dagger e^{i\delta t}]\rho\mathrm{d}W. $$ The equation is written in the interaction picture with respect to the qubit Hamiltonian $H_q = \frac{\omega}{2}\sigma_z$ and we use a detuning in the local oscillator to measure the signal appearing in the sideband of the output light.

Starting from this equation and choosing $\Delta = \omega = -\delta$, the Gaussian adiabatic elimination method [1] gives the effective dynamics $$ \mathrm{d}\rho_q = \frac{4g^2}{\kappa}\{(\bar{n}+1)\mathcal{D}[\sigma_-]+\bar{n}\mathcal{D}[\sigma_+]\}\rho_q\mathrm{d}t+ \frac{2g}{\sqrt{\kappa(2\bar{n}+1)}}\mathcal{H}[(\bar{n}+1)\sigma_-e^{-i(\phi+\pi/2)}-\bar{n}\sigma_+e^{i(\phi+\pi/2)}]\rho_q\mathrm{d}W. $$ With density operator expansion [2] for the measurement and projection operator method [3] for the deterministic part, we get $$ \mathrm{d}\rho_q = \frac{4g^2}{\kappa}\{(\bar{n}+1)\mathcal{D}[\sigma_-]+\bar{n}\mathcal{D}[\sigma_+]\}\rho_q\mathrm{d}t+ \frac{2g}{\sqrt{\kappa}}\mathcal{H}[\sigma_-e^{-i(\phi+\pi/2)}]\rho_q\mathrm{d}W. $$ See Ref. [1] for details on the derivation of these equations.

Since the full equation has time-dependent measurement operator and interaction Hamiltonian, we move to the rotating frame of the cavity field to obtain the time-independent equation $$ \mathrm{d}\rho = -ig[a\sigma_+ e^{i\phi}+a^\dagger\sigma_- e^{-i\phi},\rho]\mathrm{d}t +\kappa(\bar{n}+1)\mathcal{D}[a]\rho\mathrm{d}t+\kappa\bar{n}\mathcal{D}[a^\dagger]\rho\mathrm{d}t +\sqrt{\frac{\kappa}{2\bar{n}+1}}\mathcal{H}[(\bar{n}+1)a-\bar{n}a^\dagger]\rho\mathrm{d}W. $$ This equation is taken as the full solution with which the effective equations will be compared. The difference between the solutions is again given by the trace distance of the corresponding density operators (see the file Dispersive qubit readout.ipynb for details).

References

[1] O. Cernotik, D. Vasilyev, and K. Hammerer, arXiv:1503.07484.

[2] A. Doherty and K. Jacobs, Physical Review A 60, 2700 (1999).

[3] C. Gardiner and P. Zoller, Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer, 2004).


In [1]:
%matplotlib inline


Vendor:  Continuum Analytics, Inc.
Package: mkl
Message: trial mode expires in 29 days

In [2]:
from qutip import qeye, destroy, sigmax, sigmay, sigmaz, sigmam, sigmap, tensor, smesolve, ket2dm, basis, thermal_dm
from matplotlib.pyplot import subplots, plot
from numpy import linspace, exp, sqrt, pi

In [3]:
def trace_dist(g, phase, kappa, nbar, Ntraj, NFock, Nstep, Nsub, tmax, rhoq0):
    #***** function trace_dist *****
    #This function calculates the trace distance between the exact dynamics 
    #and the approximate solutions (Gaussian adiabatic elimination, density 
    #operator expansion) for the Jaynes-Cummings interaction.
    #Input parameters:
    #      g: qubit-cavity coupling,
    #      phase: interaction phase,
    #      kappa: intrinsic measurement (and decay) rate,
    #      nbar: thermal occupation,
    #      Ntraj: number of generated trajectories,
    #      NFock: Fock number cutoff,
    #      Nstep: number of steps for evolution,
    #      Nsub: number of substeps,
    #      tmax: maximum time,
    #      rhoq0: initial qubit state.
    #Returns:
    #      tlist: list of times,
    #      D_gae: trace distance for Gaussian adiabatic elimination,
    #      D_doe: trace distance for density operator expansion.

    #***** System operators *****
    tlist = linspace(0, tmax, Nstep)
    Id2 = qeye(2)
    IdN = qeye(NFock)
    a = tensor(Id2, destroy(NFock))
    x = (a+a.dag())/sqrt(2.)
    p = -1j*(a-a.dag())/sqrt(2.)
    sx = tensor(sigmax(), IdN)
    sy = tensor(sigmay(), IdN)
    sz = tensor(sigmaz(), IdN)
    sm = tensor(sigmam(), IdN)
    sp = tensor(sigmap(), IdN)
    
    #***** Full dynamics *****
    H_full = g*(a*sp*exp(1j*phase)+a.dag()*sm*exp(-1j*phase))
    c_op_full = [sqrt(2*kappa*nbar*(nbar+1)/(2*nbar+1))*x]
    sc_op_full = [sqrt(kappa/(2*nbar+1))*((nbar+1)*a-nbar*a.dag())]
    e_op_full = [sx, sy, sz]
    rho0 = tensor(rhoq0, thermal_dm(NFock, nbar))

    #***** Gaussian adiabatic elimination *****
    H_gae = 0*Id2
    c_op_gae = [(2*g/kappa)*sqrt(kappa*nbar*(nbar+1)/(2*nbar+1))*(sigmam()*exp(-1j*(phi+pi/2))+sigmap()*exp(1j*(phi+pi/2)))]
    sc_op_gae = [(2*g/sqrt(kappa*(2*nbar+1)))*((nbar+1)*exp(-1j*(phi+pi/2))*sigmam()-nbar*exp(1j*(phi+pi/2))*sigmap())]
    e_op_gae = [sigmax(), sigmay(), sigmaz()]

    #***** Density operator expansion *****
    H_doe = H_gae
    c_op_doe = [(2*g/sqrt(kappa))*nbar*sigmam(), (2*g/sqrt(kappa))*nbar*sigmap()]
    sc_op_doe = [(2*g/sqrt(kappa))*exp(-1j*(phi+pi/2))*sigmam()]
    e_op_doe = e_op_gae

    #***** Generating trajectories *****
    D_gae = 0
    D_doe = 0
    
    for i in xrange(0,Ntraj):
        print i
        sol_full = smesolve(H_full, rho0, tlist, c_op_full, sc_op_full, e_op_full, nsubsteps=Nsub, method='homodyne', 
                            solver='fast-milstein', store_measurement=True)
        sol_gae = smesolve(H_gae, rhoq0, tlist, c_op_gae, sc_op_gae, e_op_gae, nsubsteps=Nsub, method='homodyne', 
                           solver='fast-milstein', noise=sol_full.noise)
        sol_doe = smesolve(H_doe, rhoq0, tlist, c_op_doe, sc_op_doe, e_op_doe, nsubsteps=Nsub, method='homodyne',
                           solver='fast-milstein', noise=sol_full.noise)
        s1x = sol_full.expect[0]
        s1y = sol_full.expect[1]
        s1z = sol_full.expect[2]
        s2x = sol_gae.expect[0]
        s2y = sol_gae.expect[1]
        s2z = sol_gae.expect[2]
        s3x = sol_doe.expect[0]
        s3y = sol_doe.expect[1]
        s3z = sol_doe.expect[2]
        
        D_gae = D_gae+sqrt((s1x-s2x)**2+(s1y-s2y)**2+(s1z-s2z)**2)/2.
        D_doe = D_doe+sqrt((s1x-s3x)**2+(s1y-s3y)**2+(s1z-s3z)**2)/2.
        
    D_gae = D_gae/Ntraj
    D_doe = D_doe/Ntraj
    
    return tlist, D_gae, D_doe

Calculating average trace distance -- ensemble averaging


In [4]:
g = 0.1
phi = 0.
kappa = 1.
nbar = 0.5
Ntraj = 20
NFock = 8
Nsteps = 1000
Nsub = 250
tmax = 50

t, D_gae, D_doe = trace_dist(g, phi, kappa, nbar, Ntraj, NFock, Nsteps, Nsub, tmax, ket2dm((basis(2,0)+basis(2,1)).unit()))


0
Total run time:  15.96s
Total run time:  10.56s
Total run time:  10.25s
1
Total run time:  14.68s
Total run time:   9.67s
Total run time:   9.58s
2
Total run time:  14.72s
Total run time:   9.69s
Total run time:   9.59s
3
Total run time:  14.85s
Total run time:   9.77s
Total run time:   9.78s
4
Total run time:  14.94s
Total run time:   9.88s
Total run time:   9.85s
5
Total run time:  15.00s
Total run time:  10.08s
Total run time:   9.88s
6
Total run time:  16.08s
Total run time:  10.11s
Total run time:  10.14s
7
Total run time:  15.98s
Total run time:  10.50s
Total run time:  10.05s
8
Total run time:  14.91s
Total run time:   9.58s
Total run time:  10.24s
9
Total run time:  14.79s
Total run time:  11.07s
Total run time:  10.86s
10
Total run time:  14.59s
Total run time:   9.66s
Total run time:   9.74s
11
Total run time:  15.29s
Total run time:   9.88s
Total run time:  11.55s
12
Total run time:  14.60s
Total run time:   9.72s
Total run time:   9.60s
13
Total run time:  14.84s
Total run time:   9.65s
Total run time:   9.73s
14
Total run time:  15.29s
Total run time:  10.74s
Total run time:  10.13s
15
Total run time:  15.58s
Total run time:  10.62s
Total run time:  10.36s
16
Total run time:  15.27s
Total run time:   9.98s
Total run time:   9.81s
17
Total run time:  14.80s
Total run time:   9.70s
Total run time:   9.62s
18
Total run time:  14.36s
Total run time:   9.62s
Total run time:   9.57s
19
Total run time:  14.45s
Total run time:   9.66s
Total run time:   9.56s

In [5]:
fig, axes = subplots(figsize=(12,8))

axes.plot(t, D_gae, label='Gaussian ad. elimination')
axes.plot(t, D_doe, label='Density op. expansion')
axes.set_xlabel('Time')
axes.set_ylabel('Trace distance')
axes.legend(loc=0)


Out[5]:
<matplotlib.legend.Legend at 0x109fccd10>

Calculating average trace distance -- time averaging

Erorr with Gaussian adiabatic elimination

In [6]:
D_gae_time = D_gae.sum()/Nsteps
print D_gae_time


0.0463270610447
Error with density operator expansion

In [7]:
D_doe_time = D_doe.sum()/Nsteps
print D_doe_time


0.137791289635